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Abstract. We present a general method for Bayesian inference of the underlying covariance struc¬ 
ture of random fields on a sphere. We employ the Bipolar Spherical Harmonic (BipoSH) representation 
of general covariance structure on the sphere. We illustrate the efficacy of the method as a principled 
approach to assess violation of statistical isotropy (SI) in the sky maps of Cosmic Microwave Back¬ 
ground (CMB) fluctuations. SI violation in observed CMB maps arise due to known physical effects 
such as Doppler boost and weak lensing; yet unknown theoretical possibilities like cosmic topology and 
subtle violations of the cosmological principle, as well as, expected observational artefacts of scanning 
the sky with a non-circular beam, masking, foreground residuals, anisotropic noise, etc. We explicitly 
demonstrate the recovery of the input SI violation signals with their full statistics in simulated CMB 
maps. Our formalism easily adapts to exploring parametric physical models with non-SI covariance, 
as we illustrate for the inference of the parameters of a Doppler boosted sky map. Our approach 
promises to provide a robust quantitative evaluation of the evidence for SI violation related anomalies 
in the CMB sky by estimating the BipoSH spectra along with their complete posterior. 
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1 Introduction 

Many astronomical measurements on sky, as well as, data from other other parts of science, such as, 
geophysical and environmental modelling on the globe, deal with analysis of random fields on a sphere. 
Advancement in observational data has opened up the possibility of posing deeper inference problems 
with sophisticated analysis tools. In this paper we present a novel way to implement Bayesian inference 
of the underlying covariance structure. We develop our analysis strongly motivated by the analysis of 
cosmic microwave background (CMB) sky maps to infer the statistical isotropy (SI) of our observed 
universe. However, the methodology is applicable to Bayesian inference for other cosmological studies, 
and more widely to any study involving a scalar random field on the sphere. 

For a statistically isotropic Gaussian random field, the two-point correlation function is rota- 
tionally invariant and hence the covariance matrix of the corresponding random spherical harmonic 
coefficients, i.e., diagonal and independent of the azimuthal multipole index m. Therefore, 

the angular power spectrum {Ci = {o,*^aim)) can alone provide a complete statistical representation 
of the field on the sphere. However, in presence of SI violation, the covariance matrix (a^^o/'m') can 
depend on m and the off-diagonal components can be nonzero. This calls for a representation of the 
random field that generalizes the angular power spectrum. Bipolar Spherical harmonic representation 
(BipoSH) was first proposed by Hajian and Souradeep [1] to analyze SI violation signals on CMB. 
Any SI violation of a random field on a sphere can be completely represented in BipoSH spectra 
(A^,^^) [2-5]. The L = 0 component of the BipoSH spectra. A™, is the angular power spectrum, Ci. 

The angular power spectrum (C;) can be estimated from the observed sky map using Bayesian 
statistics. Earlier Maximum likelihood methods [6-10] were used for optimal estimation of the CMB 
power spectrum. A straight forward minimisation of the likelihood is computationally prohibitive as 
it require 0{Np^^) computational complexity, where Npi^ is the number of pixels present in the map. 
Last maximum likelihood methods on CMB dataset were used in some parts of the Boomerang and 
Maxima datasets. Since Npix increased in the modern experiments, maximum likelihood methods 
are not a viable option and pseudo-C; methods [11-13] are used since then. We can explore the full 
posterior distribution with Monte Carlo sampling. However, this approach relies of the availability 
of an efficient sampling method. At low multipoles where the posterior of the Ci are highly non- 
Gaussian due to the small number of degrees of freedom, Gibbs sampling [10] is now standard. Other 
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sampling methods, used for sampling the posterior of the CMB power spectrum, involves methods 
like Hamiltonian Monte Carlo [14] etc. 

For a SI violated Gaussian random field on a sphere, calculating the only the angular power 
spectrum by minimising the likelihood with respect to C; is not sufficient. Although, when the 
SI violation is small in comparison to the Cj, the posterior distribution of the C/ may not alter 
significantly, it is still important to minimise the likelihood with respect to the off-diagonal components 
of the covariance matrix (represented in terms of BipoSH spectra, which not only guarantees 

complete posterior statistics of Ci but also provides a complete posterior statistics of BipoSH spectra. 
We present a novel method for estimating the posterior distribution of BipoSH spectra by sampling 
the full covariance matrix. The Bayesian method for finding the posterior distribution of BipoSH 
spectra allows a complete and reliable statistical inference of the presence of SI violation over the 
sphere. For sampling the posterior, we choose the Hamiltonian Monte Carlo [14-16], a Monte Carlo 
method that uses the classical Hamiltonian mechanics for calculating the posterior distribution and 
is capable of sampling the posterior distribution faster than other conventional MCMC methods in 
such high dimensional problems. 

There has been a phenomenal improvement in CMB observations in the past few years. After 
recent data release of WMAP-9 and Planck, the observations in CMB have achieved the level of 
reliability where every single significant departure from the standard model cosmology may present a 
challenge to our understanding of the universe. The standard homogeneous and isotropic cosmological 
model dictates that the temperature fluctuation in CMB sky respect Statistical isotropy (SI). However, 
recent data from Planck hint at some putative signals of SI violation in the CMB sky [17-19]. There 
are several theoretical models proposed in literature that lead to SI violation in the CMB such as the 
non-trivial cosmic topology [6, 20-22], anisotropic cosmology, etc. [23, 24]. Weak leasing of CMB and 
Doppler boost are among the known effects that lead to SI violation in the observed CMB [25, 26]. 
Artefacts of observational reality such as anisotropic noise, foreground residuals from inadequate 
cleaning, effect of non-circular beam response function, etc. [2-4] may also yield similar SI violation. 
Therefore, it is necessary to extract the SI violation signal from observational data accounting for 
any false signal of observational systematics, and then carefully statistically assess its significance vis 
a vis peculiarity of a particular realisation of SI covariance. We apply our method of the Bayesian 
inference of random fields on sphere to analyse simulated CMB maps and recover input isotropy 
violation signals in the simulated CMB sky. Our method is shown to successfully recover the input 
signal properly up to high multipoles. Another important fact is that our Bayesian approach produces 
inferences that are data-dependent; i.e. they allow testing the relative probability of different models 
given the realized data set. This is in contrast to frequentist approaches that compare the data set 
to an ensemble of alternate data sets generated from the null hypothesis. 

The paper is organised as follows. In section 2, we provide a brief review of the BipoSH formalism 
and present the likelihood for the full covariance matrix. Section 3 discusses the details of the equation 
of motion for the HMC method. In section 4, we provide details of our handling of the computational 
issues. The analysis of simulated CMB sky maps and results are included in the section 5 of the 
paper. Here, first we do our analysis on statistically isotropic CMB maps and show that recovered 
BipoSH spectra are consistent with 0. Then we show our analysis on SI violated sky maps originated 
from WMAP beam and scan pattern, dipolar modulation and the Doppler boost signal originated 
from the motion of our galaxy respectively. In this section, we also develop and method to explicitly 
recover the Doppler boost signal by directly sampling the posterior. The final section 6 is devoted to 
discussions and conclusions. 

2 A brief review of BipoSH representation 

Consider measurements that are a linear transform of a random signal and an additive instrumental 
noise, all being fields defined on a sphere. We pose the general problem of Bayesian inference of the 
underlying covariance structure of the signal field, given the knowledge of the noise covariance and 
the linear transformation relating the measured data to the signal. 
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With no loss of any generality of the problem, we specifically consider the measurement of 
CMB temperature fluctuations in the sky. The observed sky map is a convolution of the real sky 
temperature with the instrumental beam with an addition of instrumental noise. Therefore, T( 7 ) the 
actual temperature signal of the CMB sky along the direction 7 is linearly related to the observed 
sky temperature, d{j), as 


d(7) = I i?(7, + 7V(7) , (2.1) 

where 5 ( 7 , 7 ') is the instrumental beam profile and is the instrumental noise. For a perfectly 

circular beam profile, 5 ( 7 , 7 ') = ^(7 ’ tOj assumed in this work, it is easy to deconvolve the effect of 
the beam after inferring the power spectra. However, if the beam is not circular symmetric then the 
effect of the beam depends on the full scan pattern of the experiment and its deconvolution may be 
non-trivial [3, 4, 27-29]. 

For data defined on a sphere, it is convenient to work in the spherical harmonic space. The CMB 
signal, T(^), then can be expanded in terms of spherical harmonics as 

00 I 

= E E , ( 2 . 2 ) 

l—O m— — l 

where Yimij) are the spherical harmonic functions and aim are the coefficients in the spherical har¬ 
monic basis. Similarly, the observed data, ^( 7 ), can also be expanded in spherical harmonics with 
coefficients, dim- 

For a perfectly statistically isotropic sky, the two point correlation function on sky can be ex¬ 
pressed in terms of the angular power spectrum, Q, alone as 


i^lm^l'm'') — SlmVm' — Cl6ll' Sjy 


(2.3) 


Here (...) denotes the ensemble average. However, when we allow for CMB to have SI violation. 
Cl does not provide a full description of the covariance matrix. A general covariance matrix can be 
expanded in the BipoSH representation as 


{aimO-l'm') — E J (2-4) 

771,m' 

where are the Clebsch Gordon coefficients and are the BipoSH spectra that provide a 

natural generalisation of the angular power spectrum. Since, span the entire space of the 

Covariance matrix, given a set of BipoSH spectra we can calculate the entire covariance matrix and 
vice versa. It is more convenient for most non-SI phenomena to use the even-parity BipoSH spectra 
Aflf^ = ^li^ that more closely generalise Ci. The biposh representation splits the 

covariance matrix into pieces that transform separately under the action of the group of rotation 
50(3). The scalar (rotationally invariant) term is the power spectrum, O/. 

The goal of this paper is to calculate the posterior distribution of Afjl^ from the observed sky 
map, i.e., 5(A)^)^|di^) or P{Simi'm'\dim)- Rather than computing this pdf directly we sample the joint 
probability distribution, P{Simrm', o,im\dim), and then marginalise over aim- The joint distribution 
can be obtained directly by using Bayes Theorem [ 6 , 8 , 30, 31] 


P{Slml'm'T O-lmldlm) — P{dim\aim)P{o-l 


1 


VWi 


: exp 


Iml' 

1 


n E! i^lm ^Im) ^Iml'm' i^l'm' O-l'm') 


Iml' 


/Tq-2 ^ ^Im^lml'-m'^^'m'jPi^lml'm') 

V I Iml'm' \ y Iml'm' / 


(2.5) 
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where and are the elements of the inverse matrix of Nimi'm' and Simi'm' respectively. 

P{Simi'm,') is the prior on SimVm'- The choice of P{Simi'm') has been studied in literature for the 
case of SI skymaps [10]. Here, for our analysis we use a flat prior on Simi'm’, i-e., P{Simi'm') = 1- 
However, other choices of prior such as Jeffreys prior in particular can also be used for the alanysis. 


3 Hamiltonian Monte Carlo Sampling of BipoSH spectra 

Conventional Monte Carlo techniques, such as Gibbs sampling or Metropolis Hastings, draw samples 
from a given probability distribution P{Simi'm',aim\dim) using a random walk. On the other hand, 
Hamiltonian Monte Carlo (HMC) technique based on the Classical Hamiltonian Mechanics relies on 
the fact that the density of a group of particles with random momenta placed in a potential will trace 
the potential given that all of them start from random momentum drawn from a normal distribution 
with mean 0 and co-variance M [14-16, 32, 33], where M is a positive definite matrix called the mass 
matrix and can be choosen independently. It is known that HMC method can sample the distribution 
more effectively even in very high dimensional space in comparison to other conventional MCMC 
methods. 

We sample the distribution P{Simi'm', OimMzm) using HMC with the parameters aim, ■ lu a 
Hamiltonian Monte Carlo algorithm we need to define a conjugate momentum and a mass correspond¬ 
ing to each of its parameters. We consider the conjugate momentum to aim and to be Pa,„ and 
Pj^LM and a corresponding mass respectively. The mass matrices are the positive definite 

quantity by their definition. The potential in the Hamiltonian is taken as — log P{Simi'm'■iO,im\dim) 
which leads the HMC sampler to sample the posterior of P{Simi'm',aim\dim) [16, 32]. Thus the 
Hamiltonian for the motion of this ensemble of particles is 

Erp _ p\lm 

E 7r^^-MP{Slml'm',aim\dlm))- (3.1) 

ZUln, ZuIaLM 

Im LMW 

Using Hamiltonian mechanics, the equations of motion for aim can be written as 

dlm= (3.2) 

and 


Pa,„ = 


dH 


= E ^rml'm' idt'm' ” ^Im') ” E ^Fjl'm'^Pm', 


(3.3) 


Similarly, the equations of motion for will be 


ll' ll' 


and 


(3.4) 


^ \ 

PAf^'’ = - g^LM = ^ ''imSiml'm'^l'm' 1 • 


(3.5) 


Using Eq.(2.4) and the orthogonality properties of the Clebsch Gordon coefficients [34] we obtain 


and 


( E ''tmSrJil'm'ai'm) = E 

/ mm' 


Oa^m In |5[ = tr (^SrJli.m) = E Ct^l'm'STr^^i.m' • 

S ll / rrY\'m ^ 


(3.6) 


(3.7) 
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Here tr(...) represents the trace of the enclosed matrix. 

HMC is performed in two steps. In the first step, values of the momentum variables are cho¬ 
sen from the Gaussian distribution of mean 0 and variance ttIx, where x G {aim, In the 

next step a Metropolis update is performed from the state {pai^,PA^^,^irm Afin to a new state 

{Pairn,P*ALf ,^irn, by integrating the equations of motion through a time interval of fixed size, 

At. The Hamiltonian is computed in this new state and the state is accepted with probability 
min(l, exp(—Ai?)), where AH is the change in the Hamiltonian between these two states. If the new 
state is accepted then similar operation is performed considering {a*m^ as the new position 

variable, otherwise the position is not updated from {aim, Af^J^). Choice of and itialm decides 
the stability of the integration process. HMC algorithm in general uses the Leapfrog integration algo¬ 
rithm due to its time reversal symmetry and almost symplectic nature. For Ci inference, the Leapfrog 
integration works well. However, when the the covariance matrix Simi'm' is non-diagonal, we find 
that the Leapfrog integration method diverges and needs very small step size for stable integration. 
Instead a fourth order Forest and Ruth integrator [35], which is a symplectic integrator that involves 
three Leapfrog steps, works better then a simple Leapfrog. Therefore, we use fourth order Forest and 
Ruth integrator to integrate the dynamical equations.. 

4 Computational implementation 

4.1 Limitations in L space 

The method, as described above, is applicable to Bayesian inference of the covariance structure that 

covers arbitrary SI violation. However, the finite resolution and sensitivity of measurements, pose 

information theoretic limitations to the subset of BipoSH spectra that can be expected to be inferred 

from a single sky map. For CMB sky maps that resolve up to a maximum multipole, Imux, the 

maximum number of modes of information available is ^ ^max- Hence, a systematic assessment of 

all BipoSH spectra (Lmax > L > 0) seeking to infer ^ L^^x^max coefficients would be information 

1 

limited beyond Lmax ~ ^max- Coarse grained BipoSH spectra recovered in multipole bins Al will 

allow a proportionate increase in Lmax ~ (^maxA/)3. In case of SI violation studies with a given 

model parameterisation of the BipoSH spectra, the number of parameters extracted in this approach 

together with the angular power spectrum would be similarly limited by some power of /max- For 

example, the SI violation due to weak lensing of CMB is governed by the multipoles (j)LM of the 

2 

lensing potential field. In this case, one can expect to recover (j)LM, in roughly ~ /max, independent 
bipolar multipole L bins. 

4.2 Overcoming the time complexity 

There are practical limitations arising from carrying out the computationally challenging analysis in 
reasonable time. The first practical issue that we face is related to inverting the covariance matrix. 
The need to invert the covariance matrix is required by Eqs.(3.3) and (3.7). A brute force inversion 
of the matrix is computationally prohibitive. The inversion of Eq.(3.3) can be done by using Gauss 
Seidel method. Given the expectation that, Afjl^, L ^ 0 A M ^ 0, coefficients are much smaller 
compared to the matrix Simi'm' is a diagonal dominated matrix, making it ideal for Gauss Seidel 
method. 

Brute force inversion of Simi'm' in Eq.(3.7) is difficult and time consuming except for small Imax 
(This is an important and interesting sub-case. Many of the anomalies at low multipoles, / can be 
studied with an /max of a few tens, where brute force inversion is quite possible). However, again 
using the fact that, in case of CMB signal the off-diagonal components of the matrix Simi'm' are 
expected to be much smaller than the diagonal components dominated by A™, we can use Taylor 
series expansion to invert the matrix. Simi'm' can be decomposed into Simi'm' = Dimi'm' + Oimi'm', 
where Dimi'm' is a diagonal matrix consisting only A™ part of the covariance matrix, i.e., Dimi'm' = 
^^i ^imim^ii'^rnm' and Oiml'm'is the rest of the part of the covariance matrix. Expanding {Simi'm') ^ 
into Taylor series up to the first order gives us {Simi'm') = {Dimi'm' + Oimi'm')~ = {Dimi'm')~ - 
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Figure 1. The angular power spectrum and the BipoSH coefficients up to L = 2 are inferred from a 
statistically isotropic sky map generated using HEALPix. We present the results of our analysis on two 
different realisations in the left and right columns. In the left (right) column, homogeneous, white, Gaussian 
random noise with cr„ = lOjiK (a„ = 20^K) has been added to the signal. In the top row, the solid bold 
green line shows the input power spectrum for generating the sky map. The thin dark gray line shows the 
quadratic estimate of the angular power spectrum of the particular realisation before adding noise and the 
dotted olive line is the same of the realisation after adding noise. The red data-points show the mean Ci 
recovered using joint Bayesian inference of the BipoSH coefficients up to L = 2. After the first 10 multipoles 
we plot the data-points in averages in multipole bins of Al = 20. The 2nd, 3rd and 4th row A^i and 

A^i -2 are plots for the BipoSH spectra for M = 0. The results for M yf 0 are similar and, hence, not plotted 
here. We can see that almost all the BipoSH spectra are consistent with zero within 1 to 2a, as it should be 
when the maps are drawn from a statistically isotropic covariance. It should be noted that the vertical scales 
differ in the different plots. 
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Figure 2. We show the posterior distribution of Ci, Af^ at selected multipoles I = 7, 70, 700 by 

plotting the number of sample points in different bins. Total number of sample points taken is 30,000. The 
posterior distribution is for SI sky map with noise level of = lO^K. The dashed red vertical lines mark 
the mean of each distribution. Blue dash-dot lines represent the input Ci and the BipoSH coefficients. 


{Dimi'Tn')~^ Oimi'm' ■ In realistic case, the L > 0 BipoSH coefficients being much 

smaller than Ci, this first order approximation works well in the examples studied here. 

Substituting the expressions for Dimi'm' and Oimi'm' into Eq.(3.7), we obtain 

5^ooln|5| = (2Z + l)M0° (4.1) 

and 


In |5| = (-1)^+'+''+! ^{2l+l){2V + l)Af,^/ (HOOH™,) . 


(4.2) 


These provide the set equations of motion for aim and in this approximation. These equations 
are applicable in case of weak isotropy violation, which is the case for SI violations in CMB signal. 


In cases where the SI violation signal is strong, i.e., 




0(1), the truncated Taylor expansion 


approximation used here may not hold. 


4.3 Stability of numerical integration 

Another computational issue is the choice of the numerical integration method and the mass matrix. 
In Hamiltonian integrators, though Leapfrog integrator is common because the integrator preserves 
the Hamiltonian in phase space (symplectic), the propagation error being huge we have to use a 
fourth order symplectic integrator, namely Forest and Ruth integrator, which performs better and 
the propagational errors are contained at a manageable level. 
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Forest-Ruth algorithm is a combination of three Leapfrog steps and consists of the following 


steps 


X = X + 9—v 

V = V + 9hF{x) 
h 

X = X + 9—v 


X = x + {1 — 26*)— 

-I- (1 — 29)hF{x) 


X = x + {1 — 26*)—w 


X = X + 9—v 
V = V + 9hF{x) 
nh 

X = X + 9—v 


(4.3) 


where h is the step size, 9 = {2 — ^/2) ^ and x represents the variable of integration, i.e., aim and 
in our case, and v represents the velocity, i.e.. 


P<^u 


Pj^LM 

and ———, respectively. F(x) represents 

p 

the acceleration, i.e., and —-—. It can be seen that the Forest Ruth integrator is a combination 

w 

of three Leapfrog integrator with step size 9h, (1 — 29)h and 9h, respectively. 

Choice of proper mass matrix is crucial for the stability of the integration method. If we can 
show that each of the Leapfrog step is stable then the entire integration process will also be stable. 
Here we derive the choice of the mass matrix that ensures that the integration is stable [14, 36]. 

For the equations of motion of aim we have 


= Xi+ 

Pi+i = Pi- e{S~^ + N~'^)x,+ i 
Xi+i = Xi_i_i + ^M-^pi+i 


(4.4) 


We ignore Ni^p^idpm' because that part being constant is anyway stable. The above equation can 
be written as 





1 eM-i (/ - {S-^ - 7V-1)) 


Xi 



_e(5-i_iV-i) 

1 

1 

1 

7 

7 

1 




(4.5) 


For the absolute stability of the integration process we need to ensure that the eigenvalues of the 
matrix are less than unity. The characteristic equation of this matrix depends on the choice of the 
covariance and the mass matrix. Therefore, if we choose M = (S*”^ — then the characteristic 

equation will be completely independent of the covariance matrix and it will be easy to always choose a 
step size that ensures the stability of the integration steps. However, in that case M is a non diagonal 
matrix that algebraically complicates the computation scheme. Noting again that S and N can be 
expected to be diagonally dominated, we choose the diagonal approximation, M = (^ -|- 7 ^)”^, to 
the ideal mass matrix. In this case the characteristic equation will be nearly independent of the choice 












of the mass matrix making it possible to always choose a step size for which the integration method 
stabilises. In our integrator, the step, h, is chosen such that both 9h and (1 — 29)h are less than the 
maximum value of e that is set as a requirement for the stability of the integration process. 

For the analysis and results presented in this paper, we have assumed the the noise covariance 
matrix to be diagonal in spherical harmonic space, i.e., Nimi'm’ = NiSuiSmm' ■ However, the choice of 
the mass matrix will also work for weakly anisotropic noise where the off-diagonal components of the 
noise covariance matrix is much smaller then the diagonal components. The case of masked/partial 
sky observations can also be addressed by considering a non-SI noise covariance matrix where, in pixel 
space, the variance of the noise at masked pixels is set to infinity. However, in this case the choice 
of the diagonal approximation of the mass matrix may not guarantee the stability of the integration 
process because the off-diagonal components of the noise covariance matrix in the spherical harmonic 
space could be comparable to the diagonal components. Therefore, for guaranteed stability the 
appropriate choice of the mass matrix would be non-diagonal and, hence, algebraically complicate 
the computational algorithm. In this paper, for simplicity, we restrict to full sky analysis and defer 
masked sky analysis to planned future work on observed CMB maps which may necessitate a more 
complex implementation with non-diagonal mass matrices. 

Next we discuss the choice of the mass matrix for the , where L, M ^ 0. In any realistic 
CMB map we can consider that <C Therefore, expanding Eq.(3.5) up to first order we get 

^ AT Zj at 

similar equation as of Eq.(3.3), except that the {S~^ + N~^) will be replaced by -j==^A=. There¬ 
fore, following similar logic as discussed above we choose the mass matrix for as 

The choice of the mass matrix for A™ is not directly obvious from the above arguments. We 
use the mass matrix of A™ as the inverse of its variance, which is consistent with L — 0 limit of the 
expression for the BipoSH coefficients mass matrix, and is also found to provide stable integration. 

5 Demonstration of method on simulated CMB sky maps 

In this section we demonstrate our Bayesian inference method on some representative examples. We 
consider a variety of simulated CMB maps, such as, statistically isotropic sky map, non-SI map with a 
non-circular beam signal as detected in WMAP-7 year data [17], non-SI map with dipole modulation 
signal and the Doppler boosted signal on SI violation as detected with Planck observations [19, 37, 38]. 
For all these cases, the method is tested at different noise levels. In presenting the results, we plot the 
BipoSH spectra i as defined for even parity BipoSH coefficients with 

even value of L + I + I'. 

5.1 Statistically Isotropic CMB map 

First we test the method on statistically isotropic sky maps. We produce SI map using HEALPix [39] 
at the resolution Nside = 512 pixelisation. Being statistically isotropic, all the BipoSH spectra, except 

are expected to statistically consistent with zero. We then add SI Gaussian random, white, noise 
with zero mean and standard deviation (T„ to the signal map realisation. We apply our algorithm for 
joint Bayesian inference of Ci and the other BipoSH coefficients up to L < 2. 

The results are shown in Fig. 1. We test the algorithm at two different noise levels, (T„ = lOfiK 
and (T„ = 20fj,K. For plots of Ci we adopt a hybrid axis which logarithmic up to multipole, I = 100 
and linear beyond. Also, we plot the individual mean values at each multipole, I < 10, and provide 
band power average in multipole bin size of Al = 20 at larger multipoles. The plots show that the 
algorithm recovers the input power spectrum perfectly up to the high multipoles at both the noise 
levels. We show the BipoSH coefficients, and Af ^_2 in the 2nd, 3rd and 4th column of 

Fig. 1. The values of these BipoSH coefficients are consistent with zero within 1 to 2cr. Though we 
do not present plots for the other M ^ 0 BipoSH coefficients (such as All_^,Af^ etc. ), we verify that 
they are also consistent with zero within 1 to 2cr. 

In Fig. 2 we show the marginalised probability distribution of Ci, Af^ for the multipoles 

I = 7, 70, 700 respectively for the analysis with cr„ = lOfxK noise level. In the absence of noise. 


i/(2i + l)(2C-H) 




j^LM 
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Figure 3. BipoSH spectra, and Afl’_ 2 , inferred from a non-SI sky map generated by scanning the 
sky using WMAP-9 beam and scan strategy. The dark grey plot shows the quantities obtained from a 
direct quadratic estimation from the original sky map and olive dashed curve shows the same quantities after 
adding noise with the realization. The green plot is the theoretical BipoSH values estimated by averaging the 
BipoSH samples from 30 simulations. Blue curve with red error-bars show the recovered value of the same 
from Bayesian inference. The left column corresponds to a noise level of an = lOfiK and the right column to 
a noise level of cr„ = 20/iF'.ONote that the vertical scales are different in the different plots. 


the inverse of the angular power Ci is known to follow the F-distribution that tends to Gaussian at 
high multipoles. Our analysis shows similar behaviour for the recovered posterior distributions. The 
theoretical distributions of the and are not known, although, some exploratory study has 

been carried out in literature [40]. 

5.2 Beam anisotropy sourced non-SI map 

WMAP-7 CMB map showed significant detection of SI violation signal in and [17]. Later it 
was revealed that the signal SI violation signals detected were caused due to the mild non-circularity 
of the beam response of the WMAP experiment coupled to its particular scan pattern [3]. Therefore 
as a non-SI case study, we take a SI HEALPix map and scan it with a WMAP-like scan strategy 
using the publicly available (mildly non-circular symmetric) WMAP W1 band beam response maps 
to generate a simulated time ordered data stream. From the time ordered data we reconstruct the 
map [2]. This procedure introduces SI violation signatures {Af^ and Af^^ 2 ) iii the scanned map [3, 4]. 
We then add white Gaussian noise of amplitude either lOfiK and 20/iiF to the maps. These noisy 
maps are then used to infer the BipoSH spectra. 

The BipoSH spectra of Af^ and recovered from the analysis are plotted in Fig. 3. The plot 
shows that the BipoSH spectra is very well recovered up to a very high 1. Recovery of Afj^ is very 
good at both the noise levels. We note that Afl^_ 2 , however, mildly deviates from the input signal at 
high multipoles, I > 800. 

5.3 Dipolar modulation non SI map 

In the second non-SI case study, we consider a dipole modulated sky map. We generate a SI skymap 
T( 7 ) using HEALPix and then multiply it with (.95-1-0.57^) where Td is a dipole modulation oriented 
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Figure 4. estimated from a non-SI skymap generated multiplying a SI map with (.QS + O.ST,;), where Td 

is a dipole modulation oriented along the z direction. The light gray and green plot show the quadratic 

estimations from the original skymap before and after adding noise. Blue curve with red error-bars shows the 
recovered value of the same from our inference. For this analysis we use (j„ = lOfiK. 


along z direction. The map, thus generated, has explicit SI violation signals captured in Then 

we add instrumental noise at a„ = lOfiK to the map. From this map we recover the signal in 
The results are shown in Fig. 4. It can be seen that the recovery is good at all the multipoles. 

5.4 Non-SI inference in the context of a physical model: Doppler boost 

Due to the motion of our galaxy with respect to the CMB rest frame the observed CMB signal 
experiences Doppler boost. Apart from boosting the monopole temperature, the Doppler boost also 
affects the CMB temperature fluctuations. Its effect on the temperature fluctuations are two fold. 
Firstly, it produces a modulation effect, that amplifies the temperature along the velocity direction 
and reduces in the opposite direction. Second effect is the relativistic aberration effect that squashes 
the anisotropy pattern on one side of the sky and stretches it on the other, effectively mixing the 
angular scale. As a result a specific form of SI violation is introduced in the CMB maps. Planck 2013 
results measured the non-SI signal associated with effect of the Doppler boost [19, 37]. In the final 
non-SI case study, we generate a non-SI sky map consistent with the signal from the Doppler boost 
along the z direction using CoNIGS [41]. As earlier, we add Gaussian white noise with amplitude 
10/rK and 20fiK and run our analysis on the noisy maps thus produced to recover the BipoSH spectra. 
The results of the recovery of the relevant BipoSH spectra is shown in Fig. 5. The figure shows 
the values of recovered from a particular realization at two different noise levels. For both noise 
levels, the recovered values match very well with those estimated from the input map. The recovered 
means are plotted as plotted as band power averages in multipole bins, Al = 50. 

In this particular known cause of SI violation, the non-SI signal in BipoSH is very simply related 
to the Doppler boost vector, [3 = vjc related to our peculiar motion with respect to the CMB rest 
frame.. Hence, as an illustration of non-SI parameter inference readily possible in our methodology, 
we carry out a Bayesian inference of the boost parameter, /3^^, where (3^^ = f l3Y*^^dAl is the 
spherical harmonic decomposition of /3 = -. 

The posterior of can be obtained by directly sampling the (3^^ from the probability dis¬ 
tribution , aim\dim) and considering that the only SI violation signal in the map originates 

from the Doppler boost. Under such assumption = (3 ^^^and all other BipoSH coefficients 
are 0. The equation of motion for (3^^ can be obtained as 
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Figure 5. The BipoSH estimated from a SI violation arising from Doppler boost along z direction. 

The non-SI signal map is generated using CoNIGS. The dark green plot is the input signal. Light green 
and grey plot show the from the original sky map before and after adding noise. Blue curve with red 

error-bars shows the recovered value of our HMC based inference. Left plot is for cr„ = lO^K and the right 
plot is for an = 2Q^K. The recovery of is accurate up to very high 1. 


PplM 


d In Pi^SlniVm '; \^lm) 


E 


din P{Si 

dpiM 


ml+lm '; m ) 

‘d^ 


(5.1) 


and 


plAI ^ 

TTl^iM 


(5.2) 


We take the mass parameter m^iM = 1. We can integrate these equations and infer the posterior of 
. In Fig. (6) we inferred from a SI violated map generated using CoNIGS [41]. The known 
Doppler boost injected corresponds to = —1.87 x 10“^ and (3^^ = —1.24 x 10“'*-I-1.18 x 10“^L We 
add isotropic Gaussian random noise with (t„ = 20pLK. For our analysis, we remove low multipoles and 
consider the sum of Eq.(5.1) from 201 to 1024 multipoles in the range, because at the low multipoles, 
signal is small but the errorbars being large, introduces unwarranted higher error in determining beta. 
Note that the recovery of the Doppler signal is expected to improve as one includes higher multipoles 
up to ^ ~ 2000 now available from Planck. Here, in this demonstrative example on a simulated map, 
we restrict to Imax = 1000, and also add a noise much higher than that on current CMB experiments 
to establish the method in a more adverse situation. The analysis demonstrates that models with 
small number of parameters can be inferred from high-resolution data. 



Figure 6. We calculate the posterior of /?*** from a realization generated from input = —1.87 x 10“® 
and /?** = -1.24 x 10“* -I- 1.18 x 10“^i, by directly sampling the likelihood. From the recovered posterior, we 
obtain the mean = —1.76 x 10“^ and /3** = —2.94 x 10“* -I- 0.84 x 10“®i The sample map has an isotropic 
Gaussian random noise an = 20/iiL. We generate 11,000 samples which are distributed in 20 bins in each of 
the plot. The injected value of /?**'* and the sample mean from our analysis are also marked by vertical lines 
in the graph. 
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6 Discussion and Conclnsion 


We design a general method to infer the underlying covariance structure of a random field measured 
on a sphere using a completely Bayesian technique. We employ the Bipolar Spherical Harmonic 
representation of the general covariance matrix underlying random fields on a sphere and outline the 
method for a fully Bayesian inference of the angular power spectrum and the BipoSH coefficients 
simultaneously from a single observed map. We use Hamiltonian Monte Carlo for sampling the 
posterior distributions of the BipoSH parametrization of the covariance. 

We demonstrate the method with application to simulated CMB sky maps motivated by the need 
for Bayesian assessment of the presence of the SI violation signal in observed CMB maps. We consider 
CMB sky maps that are statistically isotropic (SI), as well as, non-SI case studies that carry signatures 
of a few known cases of SI violation. We test our method at different noise levels. Our method recovers 
the angular power spectra Ci up to high multipoles, I, even in presence of large noise. The recovery of 
the BipoSH spectra is also at good fidelity at all the multipoles up to Z ~ 1000. In case of the non-SI 
Doppler boost CMB maps, we also carry out a direct Bayesian inference of the posterior distribution 
of the the governing boost paramter, /3. The recovery of the Doppler signal is expected to improve as 
one includes higher multipoles up to Z ^ 2000 now available from Planck and the appropriate noise 
level that is much lower than employed here. We expect that an application of our method to recently 
released exquisite CMB sky maps from Planck will provide reliable assessment of SI violation and 
serve as a valuable tool to assess the level of evidence for or agains the candidate anomalies. Most 
assessments of non-SI signals have to assume the angular power spectrum. Hence, it is also important 
to note that here we jointly infer the angular power spectrum together with BipoSH spectra which 
makes it specifically sensitive to non-SI phenomena that affect both, a known example of which 
is non-trivial topology of the Universe [42, 43]. Besides, hunting for cosmic signal of SI violation, 
the method provides an excellent diagnostic of non-SI residuals originating from the handling of 
unavoidable observational tasks, such as, removal of foreground emission, noise inhomogeneity and 
non symmetric beam response functions. The method should be readily extendable beyond scalar 
fields to CMB polarisation maps that are already available on the ‘full’ sky and also to shear field 
maps of future weak lensing observations over large fractions of the sky. 

The general, principled, approach to a Bayesian inference of the covariance structure in a random 
field on a sphere presented here has broad potential for application to other many aspects of cosmology 
and astronomy, as well as, more distant areas of research like geosciences and climate modelling. 
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